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ABSTRACT 

A method from molecular dynamics simulation is developed for determining local elastic 
constants of an epoxy/nanotube composite. The local values of Cn, C 33 , K 12 , and K 13 elastic 
constants are calculated for an epoxy/nanotube composite as a function of radial distance from 
the nanotube. While the results possess a significant amount of statistical uncertainty resulting 
from both the numerical analysis and the molecular fluctuations during the simulation, the 
following observations can be made. If the size of the region around the nanotube is increased 
from shells of lA to 6 A in thickness, then the scatter in the data reduces enough to observe 
trends. All the elastic constants determined are at a minimum 20A from the center of the 
nanotube. The Cn, C 33 , and K 12 follow similar trends as a function of radial distance from the 
nanotube. The K 13 decreases greater distances from the nanotube and becomes negative which 
may be a symptom of the statistical averaging. 

INTRODUCTION 

The inclusion of carbon nanotubes in structural composites for aerospace applications is 
limited to date, not only by production and processing issues, but also by their lower than 
expected mechanical perfonnance in polymeric materials. In a composite, the weak interaction 
between the nanotube and the polymer may be improved by chemical modifications 
(functionalization) to the carbon nanotube that improve the interface with the polymer. One way 
to probe the interfacial region is to examine the changes in elastic constants as a function of 
distance from the carbon nanotube. Elastic constants of polymer nanocomposites can then be 
determined for representative volume elements (RVE) of the molecular structure by using 
molecular dynamics simulation to calculate configurational energy changes under mechanical 
deformation. It is anticipated that the changes observed in the elastic constants are representative 
of the changes in molecular structure of the polymer in the vicinity of the nanotube, and 
therefore a measure of the size of the interfacial region. To understand how the nano tube affects 
the elastic constants of the polymer at the polymer nanotube interface, a broad range of spatial 
definition is required. 

Several theoretical studies have probed the ‘interfacial region’ of the polymer 
surrounding the nanotubes to elucidate the effect of the polymer-nanotube interaction on the 
mechanical behavior. Molecular dynamics (MD) simulations of nanotube pull-out and frictional 
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behavior of nanotubes in polyethylene have demonstrated the slippage of pristine nanotubes in 
polymer, and the increased resistance to motion when axial load is applied to the nanotube for 
nanotubes covalently bonded into the polymer [1, 2]. MD simulations have also been used to 
calculate the stress-strain behavior of nano tubes in polymer matrices [3]. The effect of nano tube 
functionalization on the elastic constants has been demonstrated by the equivalent continuum 
method indicating increases in the transverse elastic constants and decreases in the longitudinal 
elastic constants [4]. 

The RVEs from MD are on the order of 10,000 atoms or a few nanometers in size. Under 
deformation, they include energetic changes resulting from the chemical interactions described 
by the molecular force field. To understand how the nanotube affects the elastic constants of the 
polymer at the polymer nano tube interface, more spatial definition is required. To describe 
interfacial mechanical properties, a combined method of micromechanics techniques and 
molecular dynamics simulation was devised to determine the effects of functionalization on a 
nanoparticle interface utilizing the radial distribution function of the polymer density away from 
the nanoparticle [5]. In this study, the interfacial elastic constants were calculated for a 
functionalized silica particle in polymer for the interface defined by the region between the silica 
particle and polymer bulk density. The elastic constants were then discerned from the 
micromechanics method by comparison with elastic constants calculated for the RVE by the 
equivalent continuum method. An alternative way to calculate the interfacial elastic constants is 
as a function of radial distance C(r) from the carbon nanotube directly from the molecular 
dynamics simulation. The C(r) has been obtained directly from the atomic structure and the force 
field of a face centered cubic Lennard- Jones structure [6]. Even for Lennard- Jones interactions, 
this method is numerically complex. A different method is preferred for polymer-nanotube 
systems which additionally require consideration of bonding and Coulombic interactions. Some 
authors have used a method based on the stress fluctuations for polymers [7, 8]. Both of these 
works demonstrate that regions of negative moduli are possible within polymer systems. 

In the present work, elastic constants are calculated as a function of radial distance C(r) 
from the carbon nanotube. These local elastic constants are determined directly from the 
energies of deformation in the MD simulation for subcomponents of the RVE. The objectives of 
the present work are first to describe the method, and then to determine the spatial dependence of 
the elastic constants of the polymer surrounding a carbon nanotube as a function of distance from 
the nanotube. The spatial dependence is determined for four elastic constants the Cn, C33, K12, 
and K13 where the 3 -direction is the direction of the nano tube. The variation observed in these 
constants is significant to understanding the extent of the interfacial region in the nanocomposite. 

MATERIAL DESCRIPTION 

The representative volume element of the 
nanocomposite to be simulated contains a (10,10) 
carbon nano tube in an epoxy matrix (Figure 1). The 
epoxy matrix is composed of diglycidyl ether of 
bisphenol A (DGEBA) (Figure 2(a)) cured with 
diamino-diphenylsulfone (DDS) (Figure 2(b)). The 
molecular RVE contains about 28000 atoms of cured 

epoxy and 880 carbon atoms in the nanotube, and is Figure 1. Molecular RVE of 
approximately 10x5x5 nm in size. The epoxy is epoxy/nanotube composite. 



constructed from one 1 : 1 network of cured 
DGEBA/ DDS. It is simulated by using the 
AMBER force field[9]. Parameters from the 

literature were added for sulfur [10, 11] and 
electrostatic charges and molecular geometries were 
determined from ab initio calculations using the 
NWCHEM program[12]. The system has periodic 
boundary conditions in all 3 directions. 

The molecular RYE is simulated with 



Figure 2. Chemical structures of 
(a) DGEBA and (b) DDS. 
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molecular dynamics simulation under NVE (constant number of atoms, volume, and total 
energy) conditions using velocity rescaling to maintain a constant temperature of 300 K. The 
system is incrementally deformed every 10000 steps by a positive or negative strain of 0.0025. 
Simulations were performed for four different deformations: uniaxial strain in the 1 -direction, 
uniaxial strain in the 3-direction, plane strain in the 1- and 2-directions and plane strain in the 1- 
and 3- directions. Each of these deformations was simulated with positive or negative strain 
increments for a total of eight simulations. In each case, the strain was referenced to the 
undeformed system. The four deformations are chosen for determining the following four 
constants: Cn, C33, K12, and K13. The simulations were carried out with DL-POLY[13], 


COMPUTATIONAL MODEL 


The local elastic constants of the molecular RVE are obtained as a function of the 
distance from the nanotube. In each simulation, the configurational energy is calculated per 
atom. The configurational energy per atom A includes: 
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The V hond is one half of the potential energy associated with the chemical bond; the other half 
belongs to the bonded atom. Similarly V anf;/e is one half of the valence atom potential energy 
which is assigned to each of the two atoms connected by the angle, and V dihedral is one half the 
dihedral potential energy which is assigned to the first and fourth atoms of the dihedral. Both the 
Lennard-Jones ( V u ) and the Coulombic potential ( V Coulomh ) energies are included for atom pairs 
within a cut-off radius of 1.3 nm. The first and fourth atoms of dihedral angles contribute to both 
Vjj and V Coulomh and are scaled by 0.5 and 0.83 respectively as in the AMBER force field [9]. 

For the purposes of the analysis, no long range corrections are included for either V LJ and 
V c 0U i 0 mb ’ even though these corrections are included in the simulation. 

The A is separated into the configurational energy per ring around the nanotube, and is 
averaged over the simulation trajectory as a function of strain. Specifically, the A per ring was 
calculated from averaging the configurations every 100 time steps from 7000 time steps to 10000 
steps after the strain was incremented. The local elastic constants are then determined from the 
resulting configurational energies per ring from each of the four defonnations as a function of 
both strain and location by the following relations: 
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All the strain increments are 0.0025. For Ki 2 s u = s 12 and for K 43 s n = s It is assumed that A 0 
is the configurational energy at zero strain. The local elastic constants for each shell are 
calculated from numerically fitting A to a quadratic function of strain. The elastic constants for 
the molecular RVEs are calculated using the total configurational energies from the same 
simulations. 


RESULTS AND DISCUSSION 

The four elastic constants detennined for the Table I. Elastic Constants for 
molecular RVE are included in Table I. The uncertainty Molecular RVE 

reported is the standard error associated with the least 
squares fit of the configurational energy to the 
quadratic function of strain. The configurational energy 
varies during the simulation because the simulation 
takes place in the NVE ensemble which conserves total 
energy but allows the configurational energy to 

fluctuate. Therefore at each increment of strain the average configurational energy also varies, 
and the relationship is not perfectly quadratic. The numerical fit averages the energy fluctuations 
at different strain levels, and the standard error in the fit reflects the degree of fluctuation. Longer 
simulations may improve the results by including more configurational averaging. 

Figure 3 includes plots of the local elastic constants as a function of radial distance from 
the center of the nanotube. Results are shown for each lA shell around the nanotube, and the 
standard error in each fit is reported by the error bars. For the epoxy, the volume used is the 
volume of the shell, and for the nanotube, the volume includes the thickness of the nanotube and 
the enclosed volume. The shells are referenced with respect to the nanotube center, and do not 
change with strain. The influence of the nanotube at about 7A is expectedly strong for both the 
C33 and K13, and it barely contributes at all to C11 and K12. Any trend in epoxy behavior as a 
function of distance from the nanotube is unclear from these plots. The uncertainties in the 
elastic constants of the epoxy are much larger than the uncertainty of the elastic constants of the 
nanotube probably because of the increased mobility of the epoxy atoms relative to the nanotube, 
and therefore greater ability to sample more local environments. The negative values reflect 
situations where the energy decreased on average relative to the unstrained configuration. They 
are more likely the result of scatter in the deviations from positive moduli rather than negative 
moduli resulting from mechanically unstable parts of the material. Other authors using different 
methods have observed negative moduli when determining local elastic constants [7, 8]. 

Figure 4 includes a recalculation of the local elastic constants taking larger intervals away 
from the nanotube. In Figure 3 the local elastic constants are reported every lA, but in Figure 4, 
they are reported every 6A. The 6A interval was chosen because the first peak in the radial 
density distribution in the unstrained configuration extends from the edge of the nanotube (about 
8 A) to 14A. Therefore, it is reasonable to think of the region of epoxy up to 14A as being that of 
the first ‘ring’ of atoms surrounding the nano tube. Beyond 14 A, the radial density distribution 







Figure 3 . Radial distribution of local elastic constants. 


does not indicate clear structure so the averaging was perfonned for every 6A beyond 8 A. This 
approach avoids some of the negative moduli which are related to the unstable mechanical 
equilibrium which, in turn, may be related to how much of the group motion of the epoxy atoms 
around the nanotube is represented. It also substantially reduces the amount of scatter in the data 
making trends in the elastic constant variation within the RVE visible. All four constants show a 
minimum at 20A. Based on the uncertainty in the error bars, this minimum may not be 
statistically significant. For Cn and C33, the epoxy follows a similar trend with the elastic 
constant remaining positive. The epoxy Cn and C33 vary by about 70 GPa on average, and the 
K12 and K13 vary by about 16 and 30 GPa respectively, on average. While the K12 also remains 
positive out to 26A, the K13 is positive near the nanotube (similar to K|?) and then becomes 
negative. Ideally, the influence of the nanotube on the surrounding epoxy would be seen as 
variation in the elastic constant radially from the nanotube. In this unfunctionalized nanotube 
composite, three of the constants are unaffected by the nano tube’s presence in the epoxy to 
within the statistical uncertainty. If the trend in the Kj 3 behavior is statistically significant, then 
the nanotube affects the magnitude of K43 beyond 26A. That the K13 is negative at 20 and 26 A 
may be a symptom of either the statistical averaging or the 6A shell thickness. Other authors 
who have looked at the distribution of elastic moduli around a nanoparticle have found that there 
are regions within a ring around the nanotube where the modulus is not constant[ 7 , 8]. Further 
consideration of the size and shape of the region of the epoxy region beyond the nanotube may 
be needed to obtain a more accurate understanding of the variation of 
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Figure 4 . Radial distribution of local elastic constants at every 6A from the nanotube. 


the elastic constant with distance from the nanotube. Also, the small size of the region which 
can be properly studied with molecular dynamics simulation may be limited given the molecular 
fluctuations due to mobility of the atoms and thermal fluctuations. 

CONCLUSIONS 

A method is developed using the configurational energy within molecular dynamics 
simulations for detennining the local elastic constants of an epoxy/nanotube composite as a 
function of radial distance from the nanotube. Local Cn, C33, K12, and K [3 elastic constants are 
determined. It is found that the results depend in part upon the amount of material included in 
the analysis and that the results have significant statistical uncertainty of which is partially 
numerical and partially related to the molecular fluctuations. Increasing the size of the region 
away from the nano tube from lA shells to 6 A shells, significantly improves the scatter in the 
data to the point where trends are visible. All the elastic constants determined are at a minimum 
20A from the center of the nanotube. The Cm. C33 and K12 follow similar trends as a function of 
radial distance from the nanotube. The K13 decreases and becomes negative as the distance 
increases from the nanotube which may be a symptom of the statistical averaging. The potential 
of this method to study the elastic behavior of material in the region of the nano tube/ epoxy 
interface rests upon being able to better resolve the degree of fluctuations within the results. 
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ABSTRACT 

The inclusion of carbon nanotubes in structural composites for aerospace 
applications is limited to date, not only by production and processing issues, but also by 
their lower than expected mechanical performance in polymeric materials. In a 
composite, the weak interaction between the nanotube and the polymer may be improved 
by chemical modifications (functionalization) to the carbon nanotube that improve the 
interface with the polymer. One way to probe the interfacial region is to examine the 
changes in elastic constants as a function of distance from the carbon nanotube. Elastic 
constants of polymer nanocomposites can then be determined for representative volume 
elements (RVE) of the molecular structure by using molecular dynamics simulation to 
calculate configurational energy changes under mechanical deformation. To understand 
how the nanotube affects the elastic constants of the polymer at the polymer nanotube 
interface, a broad range of spatial definition is required. In the present work, the elastic 
constants are calculated as a function of radial distance C(r) from the carbon nanotube. 
These interfacial elastic constants are detennined directly from the energies of 
deformation in the MD simulation for subcomponents of the RVE. The objectives of the 
present work are first to describe the proposed method, and then to detennine the spatial 
dependence of the transversely isotropic elastic constants of the polymer surrounding a 
carbon nanotube as a function of distance from the nanotube. The method is applied to 
epoxy/nanotube systems in which the nanotube may be pristine or functionalized. 



